home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / obsolete / regression.pro < prev    next >
Text File  |  1997-07-08  |  6KB  |  242 lines
  1. ; $Id: regression.pro,v 1.2 1997/01/15 04:02:19 ali Exp $
  2. ;
  3. ;  Copyright (c) 1991-1997, Research Systems Inc.  All rights
  4. ;  reserved. Unauthorized reproduction prohibited.
  5.  
  6.  
  7. pro printoutr,TName,BName,SST,SSE,R,C,unit
  8.  
  9. if N_Elements( TName) EQ 0 THEN TName='Regression'
  10. if N_Elements( BName) EQ 0 THEN BName='Residual'
  11. printf,unit, " "
  12. printf,unit,  " "
  13. printf,unit,'                                ANALYSIS OF VARIANCE
  14. printf,unit," "
  15. printf,unit,'    SOURCE        SUM OF SQUARES      DF       MEAN SQUARE       F       p'
  16.       printf,unit,'*******************************************************************************'
  17.  
  18. MSSE = SSE/(R-C-1)
  19. DFE = R-C-1
  20.  
  21. MT = SST/(C)
  22.  
  23.  if MSSE NE 0 THEN FT = MT/MSSE else FT = 1.0e30
  24. printf,unit,               $
  25. Format='(A10,7X,G15.7,3X,I5,3X,G15.7,1X,G11.5,3X,F6.4)',    $
  26.        TName,SST,C,MT,FT, 1-F_Test1(FT,C,DFE)
  27.  
  28. printf,unit,Format='(A10,7X,G15.4,3X,I5,3X,G15.4)',   $
  29.  BName,SSe,DFe,MSSe
  30.  
  31. RETURN
  32. END 
  33.  
  34.  
  35.  
  36.  
  37. pro regression, X1, Y1, W1, A0, COEFF, Resid, YFit, sigma, $
  38.                 FTest, R, RMul, ChiSqr, VarNames = VarName,$
  39.                    ListName = LN,NoPrint=NP, Missing =M,   $
  40.                                            Unit = unit
  41. ;+
  42. ; NAME:
  43. ;    REGRESSION
  44. ;
  45. ; PURPOSE:
  46. ;    To augment and display the output of the library function Regress.
  47. ;    Additional output includes an anova table to test the hypothesis
  48. ;    that all coefficients are zero.A regression table is printed to the
  49. ;    screen or user specified file displaying the Coefficients, their
  50. ;    standard deviations, and  the T statistic to test for each coefficient
  51. ;    the hypothesis that it is zero.
  52. ;
  53. ; CATEGORY:
  54. ;    Statistics.
  55. ;
  56. ; CALLING SEQUENCE:
  57. ;   REGRESSION, X, Y, W,A0, COEFF, Resid, YFit, sigma, FTest, R, RMul, ChiSqr
  58. ;
  59. ; INPUTS:
  60. ;    X:    Array of independent variable data.  X must be dimensioned
  61. ;        (NTerms, NPoints) where there are Nterms coefficients to be
  62. ;        found and NPoints of sample data.  Y = column vector of NPoints
  63. ;        dependent variable values.
  64. ;
  65. ; OPTIONAL INPUTS:
  66. ;    W:    vector of weights for each equation.  For instrumental 
  67. ;        weighing, set w(i) = 1/Variance(Y(i)), for statistical 
  68. ;        weighting, w(i) = 1./Y(i) 
  69. ;
  70. ; KEYWORDS:
  71. ;     VARNAMES:    A vector of names for the independent variables to be used
  72. ;        in the output.
  73. ;                        
  74. ;      NOPRINT:    A flag to suppress printing out regression statistics.
  75. ;        The default is to print.
  76. ;       
  77. ;    LIST_NAME:    Name of output file.  Default is to the screen.
  78. ;
  79. ;      MISSING:    Missing data value.  If undefined, assume no missing data.
  80. ;        Listwise handling of missing data.
  81. ;
  82. ; OUTPUTS:
  83. ;    Anova table and regression summary written to the screen or to user
  84. ;    specified file.
  85. ;
  86. ; OPTIONAL OUTPUT PARAMETERS:
  87. ;    A0:    Constant term.
  88. ;
  89. ;    Coeff:    Vector of coefficients. Coeff has NTerm elements.
  90. ;
  91. ;    Resid:    Vector of residuals - i.e. Y - YFit.
  92. ;
  93. ;     Yfit:    Array of calculated values of Y, Npoints 
  94. ;
  95. ;    Sigma:    Vector of standard deviations for coefficients.
  96. ;
  97. ;    Ftest:    Value of F for test of fit.
  98. ;
  99. ;    Rmul:    Multiple linear correlation coefficient.
  100. ;
  101. ;    R:    Vector of linear correlation coefficient.
  102. ;
  103. ;    ChiSqr:    Reduced weighted chi square for fit.
  104. ;
  105. ; COMMON BLOCKS:
  106. ;    None.
  107. ;
  108. ; SIDE EFFECTS:
  109. ;    None.
  110. ;
  111. ; RESTRICTIONS:
  112. ;    None.
  113. ;
  114. ; PROCEDURE: 
  115. ;    See documentation for REGRESS in the User's Library.
  116. ;-
  117. ; On_Error,2
  118.  X = X1
  119.  Y = Y1
  120.  SY = size(Y)
  121.  
  122.  if N_Elements(unit) EQ 0 THEN            $     
  123.    if  N_Elements(LN)  THEN openw,unit,/get,ln else unit = -1
  124.  
  125.  X = float(X)  
  126.   
  127.  if (SY(0) lt 2) THEN         BEGIN
  128.      B = Y
  129.  ENDIF ELSE B = transpose(Y)  
  130.  
  131.  SA = size(X)
  132.  SY = Size(B)
  133.  
  134.  IF ((SA(0) eq 2) AND (SY(1) NE SA(2))) OR  $
  135.     ((SA(0) eq 1) and (Sy(1) NE SA(1))) THEN  BEGIN
  136.   printf,unit, 'regression - Incompatible arrays.'
  137.   goto,done
  138.  ENDIF
  139.  
  140.  
  141.  if N_Elements(M) THEN BEGIN
  142.  
  143.     if SA(0) eq 1 THEN BEGIN 
  144.  
  145.      here = where(X ne M, count)
  146.      if count ne 0 THEN BEGIN
  147.         X = X(here)
  148.         Y = Y(here)
  149.      ENDIF
  150.  
  151.      here = where(Y ne M, count)
  152.      if count ne 0 THEN BEGIN
  153.         X= X(here)
  154.         Y = Y(here)
  155.      ENDIF
  156.  
  157.   ENDIF  ELSE BEGIN
  158.     X = listwise([x,Y],M,rownum,rows,here)          ; removes 
  159.                                     ;cases with missing data
  160.     X = X(0:SA(1)-1,*)
  161.     B = B(here)
  162.  ENDELSE
  163.     SA = size(X)
  164.     SY = Size(B)
  165.  ENDIF
  166.  
  167.   if sa(0) eq 2 THEN siza = SA(2) else siza= SA(1)
  168.  
  169.  
  170.  if (siza LT 2) THEN BEGIN
  171.     printf,unit,"regression- need more than 1 observation
  172.     goto, DONE
  173.  ENDIF
  174.  
  175.  
  176.  
  177.  if SA(0) EQ 1 THEN  BEGIN
  178.     if N_elements(W1) ne 0 THEN W = W1
  179.     coeff=SimpRegress(X,B,W,YFit,A0,sigma,FTest,R,RMul,  $
  180.                         Chisqr) 
  181.  
  182.  ENDIF ELSE  BEGIN
  183.   if N_Elements(W1) EQ 0 THEN W=fltarr(siza)+1 else W = w1
  184.   coeff = Regress1(X,B,W,Yfit,A0,sigma,FTest,R,RMul,Chisqr   )
  185.   if N_ELEMENTS(coeff) eq 1 then  $
  186.      if coeff eq 1.e+30  THEN BEGIN
  187.         printf, unit, "Regression---Halting, correlation matrix is singular"
  188.         goto,done
  189.   ENDIF
  190.      
  191.   sigma = sigma * sqrt(Chisqr)
  192.  ENDELSE
  193.  
  194.  
  195.  if N_Elements(NP) EQ 0 THEN BEGIN
  196.   printf,unit,'Regression: Correlation Coeff =',RMul 
  197.  
  198.  
  199.  
  200.   if SA(0) eq 1 THEN VARNUM = 1 else VarNum = SA(1)
  201.   if SA(0) eq 1 then points = SA(1) else points = SA(2)
  202.  
  203.    NC=Size(VarName)
  204.    if (NC(1) NE 0 ) THEN BEGIN        ; set dependent
  205.                                       ; variable names
  206.      if  (NC(1) LT VarNum) THEN BEGIN
  207.        printf,unit,'regression-missing Variable names'
  208.        I=Indgen(varNum)
  209.        VarName=[VarName,'Var'+STRTRIM(I(NC(1):VarNum-1),2)]
  210.      ENDIF
  211.    ENDIF ELSE VarName='Var'+ STRTrim(INDgen(VarNum),2)
  212.  
  213.  
  214.    printf,unit," "
  215.    printf,unit,"  VARIABLE       COEFFICIENT            STDDEV                  T       P  "
  216.    printf,unit,format='(A10,3X,G15.7)','Constant',A0
  217.  
  218.    For i =0L,VarNum-1 DO             $
  219.  
  220.  
  221.       printf,unit,format=     $
  222.                '(A10,3X,G15.7,3x,G15.7,4x,G15.7,3x,F6.4 )',  $
  223.              VarName(i), coeff(i), sigma(i),     $
  224.              coeff(i)/sigma(i), $
  225.              2.0*(1-student1_t(abs(coeff(i)/sigma(i)),       $
  226.              points-VarNum-1))
  227.  
  228.    
  229.    Resid = B - YFit
  230.    if SA(1) NE 1 THEN ChiSqr = Total(Resid^2 * W)
  231.    Mean = Total(B)/points
  232.    SST = Total(W*(B - Mean)^2) - ChiSqr
  233.    printoutr,TName,BName,SSt,ChiSqr,Points,VarNum,unit
  234.    ChiSqr = ChiSqr/(Points - VarNum -1) 
  235.    FTEST =  SST/VarNum/ChiSqr
  236.  ENDIF
  237.  
  238.  DONE:
  239.  if (N_Elements(LN) NE 0) THEN Free_LUN,unit
  240.  return
  241.  end  
  242.